
%\documentclass[paper,twocolumn]{geophysics}
\documentclass[manuscript,revised]{geophysics}

\usepackage{multirow}

\newcommand \txt \textrm
\newcommand \vct \mathbf
\newcommand \mtx \tensor
\newcommand \ds \displaystyle

\renewcommand{\figdir}{Fig}

\begin{document}

\title{Frequency cooperative high-resolution Radon transform}

\address{
\footnotemark[1]Department of Automation,Tsinghua University \\
State Key Laboratory of Intelligent Technology and Systems\\
Tsinghua National Laboratory for Information Science and Technology \\
}


\author{Zhonghuan Chen and Wenkai Lu\footnotemark[1]}
\maketitle

\footer{Frequency cooperative high-resolution Radon transform}
\lefthead{Chen \& Lu}
\righthead{FCHRT}

\ms{FCHRT}

\begin{abstract}
In this paper, the conventional frequency-domain high-resolution Radon transform (FHRT) 
is improved in both transform performance and computation time,
by the cooperation among different frequency components.
The proposed frequency cooperative high-resolution Radon transform (FCHRT) 
retrieves the most reliable components at first,
and then uses its converged information to help other components.
Removing the iterations, FCHRT reduces the computation time greatly.
With the frequency cooperation, it protects the effective bandwidth and improves the resolution.
Both numerical and field examples are given to demonstrate these advantages.
\end{abstract}



\section{Introduction}


Radon transform (RT) has been widely applied in seismic data processing to 
separate wave fields that have different kinematic characters.
Traditionally, computation time and resolution are two of the most prevalent problems 
in field applications \cite[]{trad:386}.
In order to obtain a high resolution, 
the Radon transform is modeled to be a sparse inverse problem \cite[]{thorson:2727},
and then solved by iterative re-weighted least square (IRWLS) methods \cite[]{scales1988fast}.
%In this strategy, there is a linear equation set to be solved in every iteration.
In this strategy, a set of linear equations is used in every iteration and 
conjugate gradient optimization \cite[]{sacchi:1477,ji:R59} 
is applied to improve the computational efficiency.

For wave fields along time invariant trajectories (such as line and parabola),
a Radon transform model can be built in both the time and frequency domain.
The high-resolution solution can then be solved by IRWLS.
The time domain Radon transform (TRT) is good at expressing wave fields on complicated 
(such as time-variant) kinematic trajectories.
For example, a common-midpoint (CMP) gather can be directly mapped by hyperbolic TRT.
In the frequency domain, 
the large scale inverse problem is converted into a series of smaller problems 
for each frequency component,
%which brings a great reduction of the computation cost.
which greatly reduces the computation time.
This feature makes the frequency domain high-resolution RT (FHRT) 
more practical and popular in field applications,
despite there being an NMO-correction and an inverse NMO operation needed 
when the RT is applied to a CMP gather.

The wave fields in the TX domain consist of harmonic waves of different frequencies.
For the time-invariant Radon transform, 
each frequency component can be retrieved separately by solving 
the following equations \cite[]{trad:386}:
\begin{equation}\label{eq:hrt:model}
(\mtx W_{f}^\txt H\mtx W_{f}+\mtx L_f^\txt H\mtx W_e^\txt H\mtx W_e\mtx L_f)\bar \vct u_f=\mtx L_f^\txt H\mtx W_e^\txt H\mtx W_e\bar \vct d_f
\end{equation}
where, $\bar \vct d_f,\bar \vct u_f$ are the $f$-Hz components of 
the TX domain data and the corresponding transform domain data, 
$\mtx L_f$ is the transform matrix.
$\mtx W_f$ and $\mtx W_e$ are diagonal matrices related to the prior probability distributions 
of $\bar \vct u_f$ and approximating error $\vct e_f=\mtx L_f\bar\vct u_f-\bar\vct d_f$ separately.
The conventional Radon transform assumes that both $\bar\vct u_f$ and $\vct e_f$ are
of normal distributions with corresponding variances $\sigma_u^2$ and $\sigma_e^2$.
This means that $\mtx W_{f}=\ds{\frac{1}{\sigma_u}}\mtx I$ and 
$\mtx W_e=\ds{\frac{1}{\sigma_e}}\mtx I$.

In order to obtain a high-resolution transform, 
the transform domain $\bar\vct u_f$ is assumed to have a sparse distribution,
such as a Laplace distribution or a Cauchy distribution \cite[]{sacchi1997reweighting}.
In these cases, the diagonal weighting matrix $\mtx W_{f}$ is constructed as 
a nonlinear function of the unknown variables $\bar\vct u_f$
\begin{equation}\label{eq:kernel}
W_{f}(x,x)=g(|\bar u_f(x)|)
\end{equation}

Equation (\ref{eq:hrt:model}) and (\ref{eq:kernel}) constitute an iterative algorithm, 
known as IRWLS, proposed by \cite{scales1988fast}.
In each iteration, equation (\ref{eq:hrt:model}) becomes a linear equation,
which can be solved by Cholesky decomposition or 
the conjugate gradients optimization method \cite[]{sacchi:1477}.

%Each frequency component is retrieved by the above iterative algorithm in FHRT.
Although FHRT is more effective than the Time domain HRT (THRT), 
the computations are still quite time consuming due to the large scale pre-stacked datum.
It is common that solving the components simultaneously in THRT may obtain a better performance
than doing that separately in FHRT \cite[]{schonewille:2565}.
The experiment discussed in \cite{cary:1999} also shows that
using a bandpass filter in FHRT can obtain a better horizontal resolution.
The above research results tell us
there are beneficial relationships among these frequency components.
%Also as we know, different components of seismic data have different signal-noise-ratio (SNR),
%which leads to different iterative reliabilities.
In this paper, these relationships are exploited,
and cooperation among different frequency components is imported into FHRT.
Both theoretical analysis and straightforward experiments demonstrate the efficiency and 
the effectiveness of our proposed frequency cooperative high-resolution RT (FCHRT).
%The cooperation is based on the dispersion free assumption, but it can also work well in weak dispersion case,
%and the weak dispersion condition is present.

%Two synthetic examples are designed to show the detail

\section{Frequency cooperation}
%\old{
%Wave field information is transformed from time-offset domain into intercept-slowness
%(maybe velocity, maximum moveout or curvature according to different applications) domain in RT.
%For a frequency component, $\vct m_f$ denotes the energy distribution over different slowness of $f$-Hz harmonic waves,
%and the regularized term $\mtx W_{f}$ can be seen as a penalty for the energy distribution in the next step.}


\subsection{Dispersion free case}

In the elastic medium, both the amplitude and the shape of the wavelet are time invariable.
Considering a record of $n$ plane waves, the wave field can be expressed as:
\begin{equation}
d(t,x)=\sum_{i=0}^{n-1}h_i(t-\tau_i-p_ix)
\end{equation}
where, $h_i(t)$ is the wavelet of the $i$-th plane wave, with slowness $p_i$.
By line integrating $d(t,x)$ along a parametric curve $l$, 
we obtain the intercept-slowness domain
\begin{equation}
u(\tau,p)=\int_l d(t,x)dl=\sum_{i=0}^{n-1}\int_xh_i[\tau-\tau_i+(p-p_i)x]dx\qquad l:t=\tau+px
\end{equation}

Applying a Fourier transform along the intercept direction, the frequency-slowness domain becomes
\begin{equation}
\bar u(f,p)=\sum_{i=0}^{n-1}H_i(f)e^{-j2\pi f\tau_i}\delta(p-p_i)
\end{equation}
where, $\delta(*)$ denotes the delta function.
There are $n$ vertical lines in the frequency-slowness domain, 
each denoting a plane wave in TX domain.
%The equation means that, for the $i$-th plane wave,
%the energy of different component should be focused at the same slowness position $p_i$,
%which forms a vertical line in the frequency-slowness domain.
For the line at slowness position $p_i$,
we obtain the following amplitude relationship between two different components $f\neq f_0$:
\begin{equation}\label{eq:proportion}
|\bar u(f,p_i)|=\frac{|H_i(f)|}{|H_i(f_0)|} |\bar u(f_0,p_i)|
\end{equation}

This means $|\bar u(f,p_i)|$ and $|\bar u(f_0,p_i)|$  are proportional to each other.
In the discrete case, vectors $\bar\vct u_f$ and $\bar\vct u_{f_0}$ are proportional in amplitude.
$\bar\vct u_f$ is always normalized before applying equation (\ref{eq:kernel}),
so the converged information of component $f_0$ can be used to construct the regularization of component $f$.
%The conventional FHRT can be accelerated by frequency cooperation as:
We can calculate the dispersion free FCHRT as follows:
\begin{enumerate}
  \item Choose the reference frequency $f_0$.
  \item Calculate the $f_0$ component from equation (\ref{eq:hrt:model}) by IRWLS,
and save the corresponding converged diagonal regularization term $\mtx W_{f_0}$.
  \item Calculate other frequency components with the fixed weighting matrix $\mtx W_{f_0}$ directly.
\end{enumerate}
%Obviously, the above frequency cooperative high-resolution RT (FCHRT) can reduce the computation cost greatly.

\subsection{Reference frequency}

The performance of all components depends on the reference frequency $f_0$.
When choosing $f_0$, attention must be paid to the following three issues:
%Choosing a reference frequency is quite significant,
%because its converged information determines the performance of other components.
%The peak frequency of the data is recommended as the reference frequency $f_0$,
%with the following three :

%The first issue is considering about the distinguishing ability of events close to each other.
%First, if $f_0$ is too small,
%the resolution may be not 
\subsubsection{Resolution}

\inputdir{resolution}
In order to find the relationship between the reference frequency and the resolution,
%Considering the $f$-Hz component of a wave field with two dispersion free plane waves, written as:
considering a wavefield made of two monochromatic plane waves (of frequency $f$), 
which can be written as:
\begin{equation}\label{eq:propagation}
    d(t,x)=\sum_{i=0,1}A_i\cos[2\pi f(t-p_ix)]
\end{equation}
where $A_i$ denotes the initial amplitude at the original point $t=0,x=0$,
and $p_i(i=0,1)$ is the slowness of the corresponding plane wave.

In the time direction, 
the harmonic wave is a mixture of two cosine signal with the same frequency $f$.
In the offset direction, 
it is another mixture of two cosine signal with spatial frequency $K_i=fp_i$.
For the latter, the spatial frequency resolution is limited by the offset range 
(or mostly the maximum offset $x_M$).
From Nyquist's sampling Theorem, in order to separate the two cosine signals,
the spatial frequency interval $\Delta K=|K_0-K_1|$ should meet:
\begin{equation}
    \Delta K\geq \frac{1}{x_M}
\end{equation}

This means that only the components with frequency in the following range 
can correctly separate the two plane waves:
\begin{equation}
    f\geq \frac{1}{x_M|p_0-p_1|}
\end{equation}

Since the moveouts of the two events are $m_i=x_Mp_i$, 
the frequency range is equal to the following inequality
(for the case of parabolic RT, $m_i=x_M^2p_i$):
\begin{equation}\label{ieq:lbd}
    f\geq \frac{1}{|m_0-m_1|}
\end{equation}


\multiplot{2}{res_tx_380,res_tx_365}
{}
{Synthetic profiles with two events close to each other:
(a) Two events with slowness difference 30ms/km;
(b) Two events with slowness difference 15ms/km.
}

\tabl{parameter}{The parameters of the events in Figure~\ref{fig:res_tx_380,res_tx_365}.}
{
\begin{center}
\begin{tabular}{|c|ccccc|}
  \hline
  \multirow{2}*{PLOT} & $f_m$ & $p_0$ & $p_1$ & $\Delta p$ & $f_{L}$\\ \cline{2-6}
  & Hz & ms/km & ms/km & ms/km & Hz \\ \hline
  (a) & 35  & 350   & 380 & 30 & 33.33 \\
  (b) & 35  & 350   & 365 & 15 & 66.67 \\
  \hline
\end{tabular}
\end{center}
}

%The above inequality can be described by the following straightforward example:
A straightforward example is designed to show how FHRT and FCHRT perform 
when the above inequality holds and not.
Two noise free synthetic profiles are given in Figure \ref{fig:res_tx_380,res_tx_365}.
Each profile has two plane waves simulated by a Ricker wavelet 
with no speed dispersion nor amplitude variation.
The event parameters are given in Table \ref{tbl:parameter},
where, $f_m$ is the center frequency of the Ricker wavelet and
$f_L$ is the lower limit of the reference frequency calculated by inequality (\ref{ieq:lbd}).

\multiplot{2}{res_fm_380_cauchy,res_fm_365_cauchy,res_tm_380_cauchy,res_tm_365_cauchy}
{}
{
Transform results of profiles in Figure \ref{fig:res_tx_380,res_tx_365} by conventional FHRT :
(a) frequency-slowness domain of Profile \ref{fig:res_tx_380};
(b) frequency-slowness domain of Profile \ref{fig:res_tx_365};
(c) intercept-slowness domain of Profile \ref{fig:res_tx_380};
(d) intercept-slowness domain of Profile \ref{fig:res_tx_365}.
}

\multiplot{2}{res_fm_380_fwcauchy,res_fm_365_fwcauchy,res_tm_380_fwcauchy,res_tm_365_fwcauchy}
{}
{
Transform results of profiles in Figure \ref{fig:res_tx_380,res_tx_365} by proposed FCHRT :
(a) frequency-slowness domain of Profile \ref{fig:res_tx_380};
(b) frequency-slowness domain of Profile \ref{fig:res_tx_365};
(c) intercept-slowness domain of Profile \ref{fig:res_tx_380};
(d) intercept-slowness domain of Profile \ref{fig:res_tx_365}.
}




Figure \ref{fig:res_fm_380_cauchy,res_fm_365_cauchy,res_tm_380_cauchy,res_tm_365_cauchy} shows the transform results 
of the two profiles by conventional FHRT.
Firstly, each frequency components is retrieved with the following Cauchy kernel chosen as the nonlinear function in equation (\ref{eq:kernel}):
\begin{equation}
    g(x)=\frac{1}{x^2+\alpha^2}
\end{equation}
The results are the frequency-slowness domain as shown in Plots \ref{fig:res_fm_380_cauchy} and \ref{fig:res_fm_365_cauchy}.
Then we can get the intercept-slowness domain by inverse Fourier transform along each collumns, 
shown in Plots \ref{fig:res_tm_380_cauchy} and \ref{fig:res_tm_365_cauchy}.

In Plot \ref{fig:res_fm_380_cauchy},
energy of all components is sparse enough,
but the components with frequency lower than 33.33Hz can not 
be focused at the real slowness positions.
This agrees with the inequality (\ref{ieq:lbd}).
The inaccurate positioning of these components brings crosstalks of the two events 
in the intercept-slowness domain.
%as shown in Plot \ref{fig:res_tm_380_cauchy}.
As the two events get closer to each other in Profile \ref{fig:res_tx_365},
the problems worsen (Plots \ref{fig:res_tm_365_cauchy} and \ref{fig:res_fm_365_cauchy}).
%Plots \ref{fig:res_fm_380_cauchy,res_fm_365_cauchy,res_tm_380_cauchy,res_tm_365_cauchy}(c,d) also show that,
It is clear in both Plot \ref{fig:res_fm_380_cauchy} and Plot \ref{fig:res_fm_365_cauchy},
the slowness resolution in the low frequency band is lower than
the one in the high frequency band.
%shown in Plots \ref{fig:res_fm_380_cauchy,res_fm_365_cauchy,res_tm_380_cauchy,res_tm_365_cauchy}(c,d).
%That is to say the reference frequency $f_0$ should be as higher as possible to obtain a higher slowness resolution.


%\plot{low_freq/fwcauchy}
%{width=\textwidth}
%{Transform results of profiles in Figure \ref{fig:res_tx_380,res_tx_365} by the proposed FCHRT with $f_0=35Hz$:
%(a) intercept-slowness domain of Profile \ref{fig:res_tx_380};
%(b) intercept-slowness domain of Profile \ref{fig:res_tx_365};
%(c) frequency-slowness domain of Profile \ref{fig:res_tx_380};
%(d) frequency-slowness domain of Profile \ref{fig:res_tx_365}.
%}


The proposed FCHRT is also applied in the two synthetic profiles.
The reference frequency $f_0$ is fixed to $f_m=35Hz$, the central frequency of the wavelet.
The transform results are given in Figure \ref{fig:res_fm_380_fwcauchy,res_fm_365_fwcauchy,res_tm_380_fwcauchy,res_tm_365_fwcauchy}.
With help from the reference component, 
the components in the low frequency band can be improved 
(Plots \ref{fig:res_fm_380_fwcauchy} and \ref{fig:res_fm_365_fwcauchy}).
This removes the crosstalks and produces higher resolutions in Plots 
\ref{fig:res_tm_380_fwcauchy} and \ref{fig:res_tm_365_fwcauchy}.
%One thing to be aware of is that,
Inequality (\ref{ieq:lbd}) holds true for the two events in profile \ref{fig:res_tx_380},
but not in profile \ref{fig:res_tx_365}.
It also results in positioning errors for the two events in Plots 
\ref{fig:res_tm_365_fwcauchy} and \ref{fig:res_fm_365_fwcauchy}.

Both inequality (\ref{ieq:lbd}) and Plot \ref{fig:res_fm_365_cauchy} show that,
in order to accurately position the two events in profile \ref{fig:res_tx_365},
the reference frequency should be greater than 66.67Hz.
%and the bigger the better.
The larger the frequency, the more accurately positioned events will be.


\subsubsection{Aliasing}

Unfortunately, there is also an upper limit for $f_0$, 
which is the second issue of interest.
Assume the offset direction is sampled with interval $\Delta x$ 
(mostly 12.5m or 3.125m in practice).
Then spatial sampling frequency is $K_s=\ds{\frac{1}{\Delta x}}$.
Applying the dual form of Nyquist's sampling Theorem along the offset direction again,
the sampling frequency $K_s$ should be at least twice signals' bandwidth:
\begin{equation}
    K_s\geq 2K_i  \qquad i=0,1
\end{equation}
or there may be alias for the $i$-th plane wave at $f$-Hz.
To avoid the alias, the reference frequency should meet the following condition:
\begin{equation}\label{ieq:aliasing}
    f\leq \frac{1}{2\Delta xp_i}  \qquad i=0,1
\end{equation}

The above anti-aliasing condition can also be seen in references \cite[]{turner1990aliasing,hugonnet2001high}.
The time delay of two neighbor traces for the $i$-th linear event is given by $s(i)=\Delta xp_i$.
For a parabolic event, 
the neighbor time delay is a function of both the space position $x_k=k\Delta x$ 
and the kinematic parameter $p_i$
\begin{equation}
s(i,k)=p_i(x_k^2-x_{k-1}^2)=p_i(2k-1)\Delta^2 x
\end{equation}

We can rewrite equation (\ref{ieq:aliasing}) into the following genernal form:
\begin{equation}\label{ieq:ubd}
    f\leq \frac{1}{2s(i,k)}
\end{equation}

This means that if the above inequality does not hold at any position, 
there will be aliasing for the corresponding event (the $i$-th event).
As the number of the above equalities that do not hold increases,
the aliasing becomes worse.

\subsubsection{Signal noise ratio}

\inputdir{antinoise}

Finally, seismic signals are always polluted by random noise.
In the frequency domain, 
components in the edge band are generally more polluted than the central band.
In this case, choosing a high SNR component as the reference component 
is beneficial to improving the antinoise ability of other components.

\plot{noise_tx}{}
{
TX profile with two events in noise.
}

\multiplot{2}{noise_fm_cauchy,noise_fm_fwcauchy,noise_tm_cauchy,noise_tm_fwcauchy}
{}
{Radon transform of Profile \ref{fig:noise_tx}:
(a) frequency-slowness domain retrieved by FHRT;
(b) frequency-slowness domain retrieved by FCHRT with $f_0$=35Hz;
(c) intercept-slowness domain by FHRT;
(d) intercept-slowness domain by FCHRT.
}

%\plot{noise/noised}
%{width=\textwidth}
%{Noise example:
%(a) TX profile with two events;
%(b) intercept-slowness domain of Profile (a) retrieved by FHRT;
%(c) intercept-slowness domain retrieved by FCHRT with $f_0$=35Hz;
%(d) frequency-slowness domain by FHRT;
%(e) frequency-slowness domain by FCHRT.
%}


%A numeric example is given to show how FCHRT improve the performance by
%the cooperation of reliable component when there is strong noise.
%A numeric example is given to show how FHRT and FCHRT perform for different frequency components in noise condition.
%A numeric example is given to show how FCHRT uses the central band to help the edge band and improves the performance.
In Plot~\ref{fig:noise_tx}, 
a profile with additive background noise of standard Gaussian distribution is shown.
There are two events convoluted by a 35Hz Ricker wavelet.
The first event has a slowness of 175ms/km and a peak value SNR of 5dB.
The second event has a slowness of 200ms/km and a peak value SNR of 3dB.
The speed dispersion and amplitude variation were again not considered here.

The intercept-slowness domains retrieved by FHRT and FCHRT 
are given in Plots \ref{fig:noise_tm_cauchy} and \ref{fig:noise_tm_fwcauchy} separately.
These plots show that the energy of the two events is focused much better by FCHRT than by FHRT.
From the corresponding frequency-slowness domains in Plots \ref{fig:noise_fm_cauchy} and 
\ref{fig:noise_fm_fwcauchy}, 
it is obvious that the badly polluted components in the edge band 
can not converge at the correct slowness positions in FHRT.
FCHRT, however, can do it well, especially for the weak event.
Also, with the help of the central band, the effective bandwidth is protected better in FCHRT.
This brings a higher intercept resolution in Plot \ref{fig:noise_tm_fwcauchy}.

%From the three issues discussed above, it is clear that the peak component is the most reliable component.


\subsection{Dispersion case}

The above cooperative strategy is based on the assumption 
that there is no dispersion in the wave fields.
Inside the earth, this assumption does not hold.
In this case, the propagative slowness is a complicated function of the frequency $f$,
and equation (\ref{eq:proportion}) does not hold.
Retrieving all of the components with the constant regularization $\mtx W_{f_0}$,
the dispersion information may be damaged.

%The regularization $\mtx W_{f_0}$ applied to component $f$ can be seen as a constraint to make the energy
%to focus at the position of $p(f_0)$, the corresponding slowness of each events at the reference frequency.
The regularization $\mtx W_{f_0}$ applied to component $f$ can be seen as a constraint
to make it correspond with the reference component.
An effective method to relieve the dispersion damage is to relax this constraint 
and expand the focusing range:
Substituting equation (\ref{eq:kernel}) by the following equation can yield this result:
\begin{equation}
W_{f_0}(x,x)=g(\hat { u}_{f_0}(x))
\end{equation}
where
\begin{equation}\label{eq:window}
\hat u_{f_0}(x)=\bar u_{f_0}(x)*w(x)
\end{equation}

The retrieved component $\bar u_{f_0}(x)$ is filtered by a window function $w(x)$.
The result is then used to construct the regularization $\mtx W_{f_0}$.
Generally, a Hanning or Gaussian window is used to do this work.
The window parameter can be chosen from the tradeoff of sparsity and dispersion.

%\section{Numeric examples}
%
%As we have analyzed above, the performance in four aspects of a RT algorithm is of interest.
%They are crosstalks of closer events in low band, aliasing in high band, noises and dispersions.
%
%
%\subsection{Crosstalks of closer events in low band}
%
%
%\subsection{Aliasing in high band}
%
%\subsection{Noise case}
%
%\subsection{Dispersion case}

\section{Example}


\subsection{Numerical example}

\inputdir{schonewille}
It has been proved that, for simple model, FHRT can get a better transform result than LS solution,
but when they are applied in complex data models, FHRT is not able to provide improved resolutions
\cite[]{cary:1999,schonewille:2565}.


\plot{sch_tx}{}
{
Synthetical profile with parabolic events.
}

\multiplot{2}{sch_fm_cauchy,sch_fm_fwcauchy,sch_tm_cauchy,sch_tm_fwcauchy}
{}
{Radon transform of synthetical profile \ref{fig:sch_tx}:
(a) frequency-slowness domain retrieved by FHRT;
(b) frequency-slowness domain retrieved by FCHRT with $f_0$=50Hz;
(c) intercept-slowness domain by FHRT;
(d) intercept-slowness domain by FCHRT.
}

Figure \ref{fig:sch_tx} shows a synthetic profile 
consists of eight parabolic events with different curvatures.
As shown in Plot \ref{fig:sch_tm_cauchy},
FHRT can not get a sparse enough result, in accordance with the result 
given by \cite{schonewille:2565}.
From the frequency-moveout domain in Plot \ref{fig:sch_fm_cauchy},
we can see that, 
the components whose frequency is lower than $50Hz$ can not separate these events.
If we use a bandpass filter recommended by \cite{cary:1999}, 
which remove the low frequency components, we can get a better horizontal resolution,
but the information of low frequency will be lost.

The proposed FCHRT is also applied in this sythetic profile.
As the differences of the maximal moveout of the eight events are all $25ms$,
from inequality (\ref{ieq:lbd}), 
only components whose frequency is greater than $40Hz$ can separate these events.
The reference frequency should be greater than $40Hz$.
A higher reference frequency may work better to separate the events,
but the data is band limited and the aliasing may become worse as shown in inequality (\ref{ieq:ubd}).
Plot \ref{fig:sch_tm_fwcauchy} shows the transform result by FCHRT 
with the reference frequency $50Hz$,
which can obtain a better resolution in comparison with FHRT.
In the frequency-moveout domain (Plot \ref{fig:sch_fm_fwcauchy}), 
most of the components at lower frequencies has been well recovered.

%In this example, 30Hz Ricker wavelet is used. 
%When there is noise in the data, 
%the reference component at $50Hz$ may be difficult to converge,
%because of the low SNR at this frequency.


\subsection{Field example}

\inputdir{ec97}

\plot{cdp5200_tx_orig}
{}
{The NMO corrected CMP gather at CDP=5200.
}

\multiplot{2}{cdp5200_tx_fhrt,cdp5200_tx_fchrt}{}
{
The CMP gathher at CDP=5200 after demultiple by FHRT (a) and FCHRT (b):
The right profile has less aliasing noise than the left one at far offset area.
}

%\plot{cdp5200_tx_fchrt}{}
%{
%The CMP gathher at CDP=5200 demultiple by FCHRT:
%Compare to the far offset area in Plot \ref{fig:cdp5200_tx_fhrt}, 
%FCHRT has a better antialiasing ability.
%}

 
The two methods are applied to remove the multiples in the CMP gather of 
a 2D marine data set after NMO correction and muting.
For FCHRT, the reference frequency is chosen as the peak frequency,
and we use the 5 points Hanning window $w(n)$ to preserve the dispersion.

The CMP gather with CDP key 5200 before the demultiple is extracted 
in Figure~\ref{fig:cdp5200_tx_orig}.
The demultiple results by conventional FHRT and FCHRT are shown 
in Plots \ref{fig:cdp5200_tx_fhrt} and \ref{fig:cdp5200_tx_fchrt} respectively.
The results are quite similar, but the computation time by FCHRT is 
about one-fifth of the time by FHRT.

It is clear that there are residual multiples in both Plots \ref{fig:cdp5200_tx_fhrt} 
and \ref{fig:cdp5200_tx_fchrt} at far offset area.
The zoomed-in views tell us these residuals are results of spatial aliasing.
Comparing the two plots, FCHRT has a slightly better anti-aliasing ability.

%It becauses FCHRT used the peak component to help high frequency components.

%For the whole data set, the sampling interval in the moveout dimension is $\Delta m=0.018s$.
%The reference frequency $f_0$ is fixed as the peak frequency of the average spectrum for each CMP profile.
%And for the profile in Plot~\ref{fig:tx_cdp5200}, $f_0=35$Hz,
%and the effective frequency band of the signal is

\plot{ec97_stk_orig}
{}
{
Stacked result without demultiple.
}

\multiplot{1}{ec97_stk_fhrt,ec97_stk_fchrt}
{}
{The same parts of stacked results after demultiple by FHRT (a) and FCHRT (b).
}

Figure \ref{fig:ec97_stk_orig} shows a part of the stacked section.
Plots in Figure \ref{fig:ec97_stk_fhrt,ec97_stk_fchrt} show the demultiple results 
of the same stacked part by the FHRT and FCHRT.
In Plot \ref{fig:cdp5200_tx_orig}, 
the event at 0.4s is the primary reflection of the ocean floor,
with its first order multiple at about 0.82s.
Both the two methods have removed most of the energy of this multiple.

%but left some residuals in the sections.
%The residuals by the proposed method in Plot~\ref{fig:stk_5000_350}(c) 
%are clearly weaker than FHRT.

There is a weak slant event marked by an ellipse in the given area.
Without performing any transform, 
this event is preserved after stacking in Plot~\ref{fig:ec97_stk_orig}.
In Plot~\ref{fig:ec97_stk_fhrt}, this event has been greatly destroyed by FHRT.
In Plot~\ref{fig:ec97_stk_fchrt}, 
it is clear that the proposed method can protect it in some degree.

In Figure \ref{fig:cdp5200_tm_fhrt} and \ref{fig:cdp5200_tm_fchrt},
the transform domains of the CMP gather at CDP=5200 describe the reasons why FCHRT can protect it.
Similary to Figure \ref{fig:noise_fm_cauchy,noise_fm_fwcauchy,noise_tm_cauchy,noise_tm_fwcauchy}, 
the transform domain obtained by FCHRT has a much higher SNR than FHRT.
Zooming in on the concerned range, we can see how the slant event is focused on the intercept-moveout domain by the two methods.
At 0.82s, both FHRT and FCHRT have an energy focus for this event.
In Plot \ref{fig:cdp5200_tm_fhrt}, produced by FHRT,
this weak event has been badly overwhelmed by noise, 
while FCHRT has not (Plot \ref{fig:cdp5200_tm_fchrt}).

\plot{cdp5200_tm_fhrt}{}
{
The intercept-moveout domain of the CMP Profile at CDP=5200 by FHRT.
}

\plot{cdp5200_tm_fchrt}{}
{
The intercept-moveout domain of the CMP Profile at CDP=5200 by FCHRT.
}


\section{conclusion and discussion}

\subsection{about the computation}

Only the reference component is solved iteratively in the proposed FCHRT,
while all other components use a fixed regularized term, so the iterative process is not needed.
The computation time may be reduced to about $\ds{\frac{1}{N}}$ of the traditional FHRT,
where $N$ is the average iterations in FHRT.
This is an important advantage for large scale seismic data analysis.

To reduce the computation time in FHRT, 
the initial regularized matrix can be taken from the previous frequency slice.
The proposed FCHRT is similar, 
but it uses the most reliable frequency component instead of the previous frequency.
FCHRT has more confidence in the regularized term, so it does not use any iterations.


\subsection{about the resolution}


FHRT can only improve the slowness resolution by sparse constraints along the slowness dimension,
while THRT can improve both the intercept and the slowness resolution.
This is due to the constraints on the whole 2D domain.
Different from these methods,
as shown in Figure \ref{fig:res_fm_380_fwcauchy,res_fm_365_fwcauchy,res_tm_380_fwcauchy,res_tm_365_fwcauchy} and 
\ref{fig:noise_fm_cauchy,noise_fm_fwcauchy,noise_tm_cauchy,noise_tm_fwcauchy},
FCHRT reduces the energy leakage of frequency components at the incorrect slowness positions and
protects the effective bandwidth at the correct slowness positions.
The former improves the slowness resolution and 
the latter improves the intercept resolution.


THRT can also be seen as a combination of the models of 
equation (\ref{eq:hrt:model}) at all frequencies.
Following the experimental results by \cite{schonewille:2565},
in combination with information of all frequencies, 
THRT can obtain a sparser transform than FHRT.
Thus, THRT improves the total performance with cooperations among different components.
The proposed FCHRT just does it more explicitly.

The spatial resolution may be decreased by the window function 
imported in equation (\ref{eq:window}).
In practise, 
a Hanning window with a short length or a Gaussian window with a small variance
is enough to protect the dispersion. 
Comparing with the benefit of reducing the convergence noise,
the small effect to the spatial resolution is acceptable.
%This has a smaller effect to spatial resolution than convergence noise.

%
%but the THRT may obtain a better solution as the the reliable frequency information may help others by accident.
%
%As the regularized matrix in FHRT is from the initial non-sparse least square (LS) RT domain for each frequency slice,
%and the regular matrix in THRT is from the whole LS RT domain,
%so the resolutions of both the two directions depend on the sparseness of initial LS solution.
%the whole LS RT domain is sparser than some frequency slice,
%which makes THRT have a more efficiency than FHRT.
%The analysis in this paper agrees with this results:
%THRT builds a model of information in the whole frequency bands and solves it simultaneously,
%while FHRT builds a model for each frequency components and solves it independently.
%They are equivalent to each other at the model-building stage for time invariant models,
%but the THRT may obtain a better solution as the the reliable frequency information may help others by accident.
%That is why THRT is more efficiency than FHRT.
%Compared with THRT, FCHRT also uses the reliable frequency information to help others,
%but does it more explicitly.



\section{ACKNOWLEDGMENTS}

This work is supported by China State Key Science and Technology Projects (2011ZX05004-003 and 2011ZX05023-005-007).


\newpage
\bibliographystyle{seg}  % style file is seg.bst
\bibliography{refs}

\end{document}
